## 

library(fixest)
library(tidyverse)
library(broom)
library(sf)

## 

rm(list= ls())

## Def helper function to tidy feols output and extract n, R2

tidy_feols_extra <- function (model) {
  
  n <- model$nobs
  m_tidy <- broom::tidy(model, conf.int = F) %>% 
    mutate(conf.low = estimate - 
             qnorm(0.975) * std.error, conf.high = estimate + 
             qnorm(0.975) * std.error)
  
  # N 
  
  out <- m_tidy %>% mutate(n = n)
  
  ## R2
  
  g <- model %>% broom::glance() %>% 
    dplyr::select(matches("squared|f_"))
  out <- out %>% mutate(rsq = g$r.squared, a_rsq = g$adj.r.squared)
  
  
}

## Define make_binary

make_binary <- function (v, method = "median", q = 0.75) {
  if (!(method %in% c("mean", "median", "quantile"))) {
    (break)("Unknown method")
  }
  if (!method == "quantile") {
    cutoff = ifelse(method == "mean", mean(v, na.rm = T), 
                    median(v, na.rm = T))
    b <- ifelse(v > cutoff, 1, 0)
  }
  else {
    cutoff = quantile(v, probs = q, na.rm = T)
    b <- ifelse(v > cutoff, 1, 0)
  }
  b
}


## Dictionary


setFixest_dict(c(treated_post = "Liberalization $\\times$ Post",
                 right_total_party = 'Total right vote',
                 left_total_party = "Total left vote",
                 left_total_party_with_fdp = "Total left vote (including FDP)",
                 gdp_pc_bin = "GDP/capita",
                 gdp_pc_log = "GDP/capita"))

#### Get federal election data ### 
## Note: In BTW, right_total_party = AfD + CDU/CSU

bt <- readRDS('data/federal_elections_clean.RDS') %>%
  mutate(year = date.id) %>% 
  mutate(left_total_party_with_fdp = spd_party +
           greens_party +
           left_party + 
           fdp_party)

## Get Shapefiles 
## this is needed for conley SEs

shp_county <- read_sf("data/vg2500_krs.shp") %>% 
  dplyr::select(RS)
centr <- shp_county %>% st_centroid() %>% st_coordinates() %>% 
  data.frame() %>% dplyr::rename(lon = 1, lat = 2) %>% 
  cbind(shp_county$RS) %>% 
  dplyr::rename(county = 3)

## Add centroids to data

bt <- bt %>% 
  left_join(centr)

## Declare main data set 
## Declare main outcomes

bt_use <- bt %>% filter(year %in% c(2013, 2017))
dvars <- c("right_total_party",
           "left_total_party")

#### Table A.4: main results for federal elections ####

m1 <- fixest::feols(.[dvars] ~ post + treated + treated_post, 
                    cluster = 'AA.district.id', data = bt_use %>% 
                      filter(state == "Nordrhein-Westfalen"))
m2 <- fixest::feols(.[dvars] ~ post + treated + treated_post, 
                    cluster = 'AA.district.id', data = bt_use %>% 
                      filter(state == "Bayern"))
m3 <- fixest::feols(.[dvars] ~ post*state + treated*state + treated_post, 
                    cluster = 'AA.district.id', data = bt_use %>% 
                      filter(state %in% c("Nordrhein-Westfalen",
                                          "Bayern")))

n_models <- 2

## Tidy these

l1 <- lapply(1:n_models, function(i) tidy_feols_extra(m1[[i]])) %>% 
  reduce(rbind)
l2 <- lapply(1:n_models, function(i) tidy_feols_extra(m2[[i]])) %>% 
  reduce(rbind)
l3 <- lapply(1:n_models, function(i) tidy_feols_extra(m3[[i]])) %>% 
  reduce(rbind)

## Get proper names from dict

d <- data.frame(getFixest_dict()) %>% 
  cbind(rownames(.)) %>% 
  dplyr::select(outcomes = 2, outcome_lab = 1)

## Prep for table output

l_out <- bind_rows(l1, l2, l3) %>% 
  filter(term == 'treated_post') %>% 
  mutate(outcomes = rep(dvars, 3)) %>% 
  left_join(d) %>% 
  dplyr::select(outcome_lab, estimate, 
                std.error, p.value, n, rsq) %>% 
  mutate(state = rep(c("NRW", "BY", "NRW+BY"), each = n_models)) %>% 
  arrange(state, outcome_lab) %>% 
  mutate(across(2:6, ~round(., 3))) %>% 
  dplyr::select(outcome_lab, state, everything())

## Kable 

library(kableExtra)

l_out %>%
  kableExtra::kable(., format = 'latex', 
                    booktabs = T,
                    linesep = "",
                    caption = 'Details on main results',
                    col.names = c('Outcome', 'Sample',
                                  'Estimate', 'SE', 'p-value', 'N',
                                  'R-squared')) %>%
  kable_styling(latex_options = c("hold_position", 'scale_down'), 
                position = "center") 

#### Table A.5: main results for state elections ####

# Bayern:
# totright2_nofdp = CDU/CSU + AfD 
# totleft2 = SPD and Greens (note that Left is part of other for BY)

# NRW:
# left_total = SPD + Greens + Left
# totright2_nofdp = CDU/CSU + AfD 

## Get Bayern data 

by <- read_rds("data/bavaria_elections_clean.rds") %>% 
  filter(year %in% c(2013, 2018)) %>%
  dplyr::select(-year) %>%
  mutate(turnout = turnout / 100) %>% 
  mutate(totright2_nofdp = CSU_2 + ifelse(post == 1, AFD_2, 0)) %>% 
  mutate(totright2_nofdp = totright2_nofdp * 100,
         totleft2 = totleft2 * 100)

## Get NRW data 

nrw <- read_rds("data/nrw_elections_clean.rds")

## Declare outcomes 

outcomes_by <- c("totright2_nofdp", "totleft2")
outcomes_nrw <- c("totright2_nofdp", "left_total")

## Estimate models

m1_nrw <- fixest::feols(.[outcomes_nrw] ~ post + treated + treated_post, 
                        cluster = 'AA.district.id', data = nrw)
m2_by <- fixest::feols(.[outcomes_by] ~ post + treated + treated_post, 
                       cluster = 'AA.district.id', data = by)

## Number of models per object 

n_models <- 2

## Tidy these

l1 <- lapply(1:n_models, function(i) tidy_feols_extra(m1_nrw[[i]])) %>% 
  reduce(rbind) %>% 
  mutate(sample = 'NRW') %>% 
  filter(term == "treated_post") %>% 
  mutate(dv = outcomes_nrw)
l2 <- lapply(1:n_models, function(i) tidy_feols_extra(m2_by[[i]])) %>% 
  reduce(rbind) %>% 
  mutate(sample = 'BY') %>% 
  filter(term == "treated_post") %>% 
  mutate(dv = outcomes_by)

## Prep for kable

ltw <- rbind(l1, l2) %>% 
  mutate(dv = recode(dv, 
                     `left_total` = 'Total left vote',
                     `totright2_nofdp` = 'Total right vote',
                     `totleft2` = 'Total left vote')) %>%  
  dplyr::select(dv, sample, estimate, std.error, p.value, n, rsq) %>% 
  arrange(sample, dv)

## ## Return this as a table 

ltw %>%
  kableExtra::kable(., format = 'latex', 
                    digits =3,
                    booktabs = T,
                    linesep = "",
                    caption = 'Details on main results',
                    col.names = c('Outcome', 'Sample',
                                  'Estimate (p.p.)', 'SE', 'p-value', 'N',
                                  'R-squared')) %>%
  kable_styling(latex_options = c("hold_position", 'scale_down'), 
                position = "center") 

#### Table A.6: Main results, adjusting standard errors for spatial dependence ####

m1_c <- fixest::feols(right_total_party ~ post + treated + treated_post, 
                      data = bt_use %>% 
                        filter(state == "Nordrhein-Westfalen"),
                      vcov = vcov_conley(lat = "lat", lon = "lon", 
                                         cutoff = 100, 
                                         distance = "spherical"))

m2_c <- fixest::feols(right_total_party ~ post + treated + treated_post, 
                      data = bt_use %>% 
                        filter(state == "Bayern"),
                      vcov = vcov_conley(lat = "lat", lon = "lon", 
                                         cutoff = 100, 
                                         distance = "spherical"))
m3_c <- fixest::feols(right_total_party ~ post*state + treated*state + treated_post, 
                      data = bt_use %>% 
                        filter(state %in% c("Nordrhein-Westfalen",
                                            "Bayern")),
                      vcov = vcov_conley(lat = "lat", lon = "lon", 
                                         cutoff = 100, 
                                         distance = "spherical"))

## Table


mlist <- list(m1, 
              m2, 
              m3, 
              m1_c, 
              m2_c, 
              m3_c)

etable(mlist, style.tex = style.tex("qje"), 
       fitstat = ~ r2 + n, tex = TRUE, keep = 'Liberalization', 
       title = "Main results, adjusting standard errors for spatial dependence", 
       placement = '!h', label = 'tab:conley', 
       digits = 4, digits.stats = 3)

#### Table A.7: Main results, controlling for time varying GDP/capita and unemployment rates #### 

## First, create binary and log version of GDP/pc

bt_use <- bt_use %>% 
  mutate(gdp_pc = as.numeric(gdp_pc)) %>% 
  group_by(year) %>% 
  mutate(gdp_pc_bin = make_binary(gdp_pc),
         gdp_pc_log = log(gdp_pc)) %>% 
  ungroup()

##

m1_gdp <- fixest::feols(right_total_party ~ post + treated + treated_post +
                          gdp_pc_bin + unem_rate, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state == "Nordrhein-Westfalen"))
m2_gdp <- fixest::feols(right_total_party ~ post + treated + treated_post+
                          gdp_pc_bin + unem_rate, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state == "Bayern"))
m3_gdp <- fixest::feols(right_total_party ~ post*state + treated*state + treated_post+
                          gdp_pc_bin + unem_rate, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state %in% c("Nordrhein-Westfalen",
                                              "Bayern")))
m4_gdp <- fixest::feols(right_total_party ~ post + treated + treated_post +
                          gdp_pc_log + unem_rate, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state == "Nordrhein-Westfalen"))
m5_gdp <- fixest::feols(right_total_party ~ post + treated + treated_post+
                          gdp_pc_log + unem_rate, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state == "Bayern"))
m6_gdp <- fixest::feols(right_total_party ~ post*state + treated*state + treated_post+
                          gdp_pc_log + unem_rate, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state %in% c("Nordrhein-Westfalen",
                                              "Bayern")))

## Table 

mlist <- list(m1_gdp, m2_gdp, m3_gdp, m4_gdp ,m5_gdp ,m6_gdp)

etable(mlist, style.tex = style.tex("qje"), 
       fitstat = ~ r2 + n, tex = TRUE, keep = c('Liberalization'),
       drop = 'state',
       title = "Main results, controlling for GDP", 
       placement = '!h', label = 'tab:unem', digits = 4, digits.stats = 3)

#### Table A.8: Main results, weighting municipalities by population ####

m1_weight <- fixest::feols(right_total_party ~ post + treated + treated_post, 
                           cluster = 'AA.district.id', data = bt_use %>% 
                             filter(state == "Nordrhein-Westfalen"),
                           weights = ~n_votes_party)
m2_weight <- fixest::feols(right_total_party ~ post + treated + treated_post, 
                           cluster = 'AA.district.id', data = bt_use %>% 
                             filter(state == "Bayern"),
                           weights = ~n_votes_party)
m3_weight <- fixest::feols(right_total_party ~ post*state + treated*state + treated_post, 
                           cluster = 'AA.district.id', data = bt_use %>% 
                             filter(state %in% c("Nordrhein-Westfalen",
                                                 "Bayern")),
                           weights = ~n_votes_party)

## Table

mlist <- list(m1, 
              m2, 
              m3, 
              m1_weight, 
              m2_weight, 
              m3_weight)

etable(mlist, style.tex = style.tex("qje"), 
       fitstat = ~ r2 + n, tex = TRUE, keep = 'Liberalization', 
       title = "Main results, weighting municipalities by population", 
       placement = '!h', label = 'tab:weight', digits = 4, digits.stats = 3)

#### Table A.9: including FDP in left vote outcome ####

m1_fdp <- fixest::feols(left_total_party_with_fdp ~ post + treated + 
                          treated_post, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state == "Nordrhein-Westfalen"))
m2_fdp <- fixest::feols(left_total_party_with_fdp ~ post + treated + 
                          treated_post, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state == "Bayern"))
m3_fdp <- fixest::feols(left_total_party_with_fdp ~ post*state + 
                          treated*state + treated_post, 
                        cluster = 'AA.district.id', data = bt_use %>% 
                          filter(state %in% c("Nordrhein-Westfalen",
                                              "Bayern")))

## Table

mlist <- list(m1_fdp, m2_fdp, m3_fdp)

etable(mlist, style.tex = style.tex("qje"), 
       fitstat = ~ r2 + n, tex = TRUE, keep = 'Liberalization', 
       title = "Main results, alternative outcome definition", 
       placement = '!h', label = 'tab:alt_fdp', 
       digits = 4, digits.stats = 3)

#### Table A.10: Effects on incumbent vote share ####

## Define incumbent outcomes

bt_use <- bt_use %>% 
  mutate(inc_state = ifelse(state == 'Bayern', cdu_csu_party, 
                            spd_party + greens_party),
         inc_joint =  ifelse(state == 'Bayern', cdu_csu_party + spd_party, 
                             cdu_csu_party + spd_party + greens_party),
         inc_fed = cdu_csu_party + spd_party)

## Estimate 

m1_inc <- feols(inc_state ~ post  + treated + treated_post,
                cluster = 'AA.district.id',
                data = subset(bt_use, state == 'Bayern'))

m2_inc <- feols(inc_state ~ post  + treated + treated_post,
                cluster = 'AA.district.id',
                data = subset(bt_use, state != 'Bayern'))

m3_inc <- feols(inc_fed ~post  + treated + treated_post,
                cluster = 'AA.district.id',
                data = subset(bt_use, state == 'Bayern'))

m4_inc <- feols(inc_fed ~ post  + treated + treated_post,
                cluster = 'AA.district.id',
                data = subset(bt_use, state != 'Bayern'))

m5_inc <- feols(inc_joint ~ post  + treated + treated_post,
                cluster = 'AA.district.id',
                data = subset(bt_use, state == 'Bayern'))

m6_inc <- feols(inc_joint ~ post  + treated + treated_post,
                cluster = 'AA.district.id',
                data = subset(bt_use, state != 'Bayern'))


mlist <- list(m1_inc, m2_inc, m3_inc, m4_inc, m5_inc, m6_inc)

etable(mlist, style.tex = style.tex("qje"), 
       fitstat = ~ r2 + n, tex = TRUE, keep = c('Liberalization'),
       title = "Effects on incumbent vote share", 
       placement = '!h', 
       label = 'tab:incumbency', 
       digits = 4, 
       digits.stats = 3)

#### Table A.15: Main results with district-specific time trends ####

## We use elections in 2009, 2013 and 2017
## Years enter in a linear fashion

bt_long <- bt %>% 
  mutate(period = as.numeric(as.factor(year)))

## Models

m1_long <- fixest::feols(right_total_party ~ treated_post +
                           AA.district.id:period| 
                           ags + period,
                         cluster = 'AA.district.id', data = bt_long %>% 
                           filter(state %in% c("Nordrhein-Westfalen",
                                               "Bayern")) %>% 
                           filter(date.id > 2005))
m2_long <- fixest::feols(right_total_party ~ treated_post +
                           AA.district.id:period| 
                           ags + period^state,
                         cluster = 'AA.district.id', data = bt_long %>% 
                           filter(state %in% c("Nordrhein-Westfalen",
                                               "Bayern")) %>% 
                           filter(date.id > 2005))

# Table

mlist <- list(m1_long, m2_long)

## ##

etable(mlist,  tex = T,
       fitstat = ~ r2 + n, keep = c('Liberalization'),
       drop = 'state',
       title = "Main results with district-specific time trends", 
       placement = '!h', label = 'tab:trend', digits = 4, digits.stats = 3 )


